Big Mac Index Geo Network Analysis

Take Home Exercise 04

This exercise aims to develop a Shiny application featuring an interactive choropleth map to analyze global trade networks in conjunction with the Big Mac Index. By integrating geographical data with economic indicators, the application will enable a visual exploration of how trade volumes and net exports correlate with price parity across different regions. The choropleth will serve as a dynamic tool to observe patterns, offering users the ability to drill down into country-specific trade connections and Index values.

The final deliverable will be a user-centric Shiny application, optimized for engagement and insight discovery. It will be complemented by a concise report summarizing the analytical narratives that emerge from the visualization. This endeavor will not only underscore the practicality of geospatial data in economic analysis but also aims to enrich the discourse on the implications of trade dynamics for market pricing mechanisms like the Big Mac Index.

1 Load Packages

For this Take Home Exercise, the focus will be on packages such as ggforce, igraph, ggraph and visNetwork.

Code

library(tidyverse)
library(sf)
library(rnaturalearth)
library(countrycode)
library(ggrepel)
library(plotly)
library(cartogram)
library(ggforce)
library(igraph)
library(dplyr)
library(png)
library(visNetwork)
library(htmlwidgets)
library(base64enc)
library(tidygraph)
library(ggplot2)
library(ggraph)
library(tmap)
library(leaflet)

2 EDA Trade on GEO

2.1 Appreciating through Choropleth

Code

choropleth<- st_read("C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/geospatial/TM_WORLD_BORDERS-0.3.shp")

Reading layer `TM_WORLD_BORDERS-0.3' from data source 
  `C:\FirGhaz\ISSS608-VAA\Take-home_Exercises\Take-Home_Ex04\data\geospatial\TM_WORLD_BORDERS-0.3.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 246 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6236
Geodetic CRS:  WGS 84

Code

choropleth_projected <- st_transform(choropleth, "+proj=robin")  
dataV <- read_csv("C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/vic_countries_with_complete_data.csv")
invalid_geoms <- st_is_valid(choropleth_projected, reason = TRUE)
print(invalid_geoms[invalid_geoms != "Valid"])

  [1] "Valid Geometry"                                             
  [2] "Valid Geometry"                                             
  [3] "Valid Geometry"                                             
  [4] "Valid Geometry"                                             
  [5] "Valid Geometry"                                             
  [6] "Valid Geometry"                                             
  [7] "Valid Geometry"                                             
  [8] "Valid Geometry"                                             
  [9] "Valid Geometry"                                             
 [10] "Valid Geometry"                                             
 [11] "Valid Geometry"                                             
 [12] "Valid Geometry"                                             
 [13] "Valid Geometry"                                             
 [14] "Valid Geometry"                                             
 [15] "Valid Geometry"                                             
 [16] "Valid Geometry"                                             
 [17] "Valid Geometry"                                             
 [18] "Valid Geometry"                                             
 [19] "Valid Geometry"                                             
 [20] "Valid Geometry"                                             
 [21] "Valid Geometry"                                             
 [22] "Valid Geometry"                                             
 [23] "Valid Geometry"                                             
 [24] "Ring Self-intersection[-5077424.46853944 7124752.10050202]" 
 [25] "Valid Geometry"                                             
 [26] "Valid Geometry"                                             
 [27] "Valid Geometry"                                             
 [28] "Valid Geometry"                                             
 [29] "Valid Geometry"                                             
 [30] "Valid Geometry"                                             
 [31] "Valid Geometry"                                             
 [32] "Valid Geometry"                                             
 [33] "Ring Self-intersection[-5608127.39693509 -5808904.89708338]"
 [34] "Valid Geometry"                                             
 [35] "Valid Geometry"                                             
 [36] "Valid Geometry"                                             
 [37] "Valid Geometry"                                             
 [38] "Valid Geometry"                                             
 [39] "Valid Geometry"                                             
 [40] "Valid Geometry"                                             
 [41] "Valid Geometry"                                             
 [42] "Valid Geometry"                                             
 [43] "Valid Geometry"                                             
 [44] "Valid Geometry"                                             
 [45] "Valid Geometry"                                             
 [46] "Valid Geometry"                                             
 [47] "Valid Geometry"                                             
 [48] "Valid Geometry"                                             
 [49] "Valid Geometry"                                             
 [50] "Valid Geometry"                                             
 [51] "Valid Geometry"                                             
 [52] "Valid Geometry"                                             
 [53] "Valid Geometry"                                             
 [54] "Valid Geometry"                                             
 [55] "Valid Geometry"                                             
 [56] "Valid Geometry"                                             
 [57] "Valid Geometry"                                             
 [58] "Valid Geometry"                                             
 [59] "Valid Geometry"                                             
 [60] "Valid Geometry"                                             
 [61] "Valid Geometry"                                             
 [62] "Valid Geometry"                                             
 [63] "Valid Geometry"                                             
 [64] "Valid Geometry"                                             
 [65] "Valid Geometry"                                             
 [66] "Valid Geometry"                                             
 [67] "Valid Geometry"                                             
 [68] "Valid Geometry"                                             
 [69] "Valid Geometry"                                             
 [70] "Valid Geometry"                                             
 [71] "Valid Geometry"                                             
 [72] "Valid Geometry"                                             
 [73] "Valid Geometry"                                             
 [74] "Valid Geometry"                                             
 [75] "Valid Geometry"                                             
 [76] "Valid Geometry"                                             
 [77] "Valid Geometry"                                             
 [78] "Valid Geometry"                                             
 [79] "Valid Geometry"                                             
 [80] "Valid Geometry"                                             
 [81] "Valid Geometry"                                             
 [82] "Valid Geometry"                                             
 [83] "Valid Geometry"                                             
 [84] "Valid Geometry"                                             
 [85] "Valid Geometry"                                             
 [86] "Valid Geometry"                                             
 [87] "Valid Geometry"                                             
 [88] "Valid Geometry"                                             
 [89] "Valid Geometry"                                             
 [90] "Valid Geometry"                                             
 [91] "Valid Geometry"                                             
 [92] "Valid Geometry"                                             
 [93] "Valid Geometry"                                             
 [94] "Valid Geometry"                                             
 [95] "Valid Geometry"                                             
 [96] "Valid Geometry"                                             
 [97] "Valid Geometry"                                             
 [98] "Valid Geometry"                                             
 [99] "Valid Geometry"                                             
[100] "Valid Geometry"                                             
[101] "Valid Geometry"                                             
[102] "Valid Geometry"                                             
[103] "Valid Geometry"                                             
[104] "Valid Geometry"                                             
[105] "Valid Geometry"                                             
[106] "Valid Geometry"                                             
[107] "Valid Geometry"                                             
[108] "Valid Geometry"                                             
[109] "Valid Geometry"                                             
[110] "Valid Geometry"                                             
[111] "Valid Geometry"                                             
[112] "Valid Geometry"                                             
[113] "Valid Geometry"                                             
[114] "Valid Geometry"                                             
[115] "Valid Geometry"                                             
[116] "Valid Geometry"                                             
[117] "Valid Geometry"                                             
[118] "Valid Geometry"                                             
[119] "Valid Geometry"                                             
[120] "Valid Geometry"                                             
[121] "Valid Geometry"                                             
[122] "Valid Geometry"                                             
[123] "Valid Geometry"                                             
[124] "Valid Geometry"                                             
[125] "Valid Geometry"                                             
[126] "Valid Geometry"                                             
[127] "Valid Geometry"                                             
[128] "Valid Geometry"                                             
[129] "Valid Geometry"                                             
[130] "Valid Geometry"                                             
[131] "Valid Geometry"                                             
[132] "Valid Geometry"                                             
[133] "Valid Geometry"                                             
[134] "Valid Geometry"                                             
[135] "Valid Geometry"                                             
[136] "Valid Geometry"                                             
[137] "Valid Geometry"                                             
[138] "Valid Geometry"                                             
[139] "Valid Geometry"                                             
[140] "Valid Geometry"                                             
[141] "Valid Geometry"                                             
[142] "Valid Geometry"                                             
[143] "Valid Geometry"                                             
[144] "Valid Geometry"                                             
[145] "Valid Geometry"                                             
[146] "Valid Geometry"                                             
[147] "Valid Geometry"                                             
[148] "Valid Geometry"                                             
[149] "Valid Geometry"                                             
[150] "Valid Geometry"                                             
[151] "Valid Geometry"                                             
[152] "Valid Geometry"                                             
[153] "Valid Geometry"                                             
[154] "Valid Geometry"                                             
[155] "Ring Self-intersection[396540.032895809 6491128.910613]"    
[156] "Valid Geometry"                                             
[157] "Valid Geometry"                                             
[158] "Valid Geometry"                                             
[159] "Valid Geometry"                                             
[160] "Valid Geometry"                                             
[161] "Valid Geometry"                                             
[162] "Valid Geometry"                                             
[163] "Valid Geometry"                                             
[164] "Valid Geometry"                                             
[165] "Valid Geometry"                                             
[166] "Valid Geometry"                                             
[167] "Valid Geometry"                                             
[168] "Valid Geometry"                                             
[169] "Valid Geometry"                                             
[170] "Valid Geometry"                                             
[171] "Valid Geometry"                                             
[172] "Valid Geometry"                                             
[173] "Valid Geometry"                                             
[174] "Valid Geometry"                                             
[175] "Ring Self-intersection[11835947.0260938 5255647.69546204]"  
[176] "Valid Geometry"                                             
[177] "Valid Geometry"                                             
[178] "Valid Geometry"                                             
[179] "Valid Geometry"                                             
[180] "Valid Geometry"                                             
[181] "Valid Geometry"                                             
[182] "Valid Geometry"                                             
[183] "Valid Geometry"                                             
[184] "Valid Geometry"                                             
[185] "Valid Geometry"                                             
[186] "Valid Geometry"                                             
[187] "Valid Geometry"                                             
[188] "Valid Geometry"                                             
[189] "Valid Geometry"                                             
[190] "Valid Geometry"                                             
[191] "Valid Geometry"                                             
[192] "Valid Geometry"                                             
[193] "Valid Geometry"                                             
[194] "Valid Geometry"                                             
[195] "Valid Geometry"                                             
[196] "Valid Geometry"                                             
[197] "Valid Geometry"                                             
[198] "Valid Geometry"                                             
[199] "Valid Geometry"                                             
[200] "Valid Geometry"                                             
[201] "Valid Geometry"                                             
[202] "Valid Geometry"                                             
[203] "Valid Geometry"                                             
[204] "Valid Geometry"                                             
[205] "Valid Geometry"                                             
[206] "Valid Geometry"                                             
[207] "Valid Geometry"                                             
[208] "Valid Geometry"                                             
[209] "Valid Geometry"                                             
[210] "Valid Geometry"                                             
[211] "Valid Geometry"                                             
[212] "Valid Geometry"                                             
[213] "Valid Geometry"                                             
[214] "Valid Geometry"                                             
[215] "Valid Geometry"                                             
[216] "Valid Geometry"                                             
[217] "Valid Geometry"                                             
[218] "Valid Geometry"                                             
[219] "Valid Geometry"                                             
[220] "Valid Geometry"                                             
[221] "Valid Geometry"                                             
[222] "Valid Geometry"                                             
[223] "Valid Geometry"                                             
[224] "Valid Geometry"                                             
[225] "Valid Geometry"                                             
[226] "Valid Geometry"                                             
[227] "Valid Geometry"                                             
[228] "Valid Geometry"                                             
[229] "Valid Geometry"                                             
[230] "Valid Geometry"                                             
[231] "Valid Geometry"                                             
[232] "Valid Geometry"                                             
[233] "Valid Geometry"                                             
[234] "Valid Geometry"                                             
[235] "Valid Geometry"                                             
[236] "Valid Geometry"                                             
[237] "Valid Geometry"                                             
[238] "Valid Geometry"                                             
[239] "Valid Geometry"                                             
[240] "Valid Geometry"                                             
[241] "Valid Geometry"                                             
[242] "Valid Geometry"                                             
[243] "Valid Geometry"                                             
[244] "Valid Geometry"                                             
[245] "Valid Geometry"                                             
[246] "Valid Geometry"                                             

Code

choropleth_projected <- st_make_valid(choropleth_projected)

Code

cleaned_dataV <- dataV %>%
  filter(!is.na(country) & !is.na(year))

map_for_V <- choropleth_projected%>%
  left_join(cleaned_dataV, by = c("ISO3" = "ISO3", "LON" = "LON", "LAT"= "LAT"))
  #select(c(country:bmi_diff))%>%
  #mutate(across(c(bmi_usdprice:bmi_diff), ~ifelse(is.na(.), 0, .)))

library(tmap)
map_for_V_2021 <- map_for_V %>%
  filter(year == 2021)
map_for_V_2021$iso2 <- tolower(map_for_V_2021$ISO2)
tmap_mode("view")  

tm <- tm_shape(map_for_V_2021) +
  tm_polygons("bmi_usdprice",  # The variable to visualize
              palette = "YlOrBr",  # Color palette
              title = "BMI (US$)",  # Legend title
              id = "country") +  
  tm_basemap("CartoDB.DarkMatterNoLabels") + # Tooltip: country name from your dataset
  tm_layout(legend.position = c("left", "bottom"))  # Adjust legend position

# Print the map
print(tm)

Looking through the lens of Continents

Code

map_for_V_2021 <- map_for_V %>%
  filter(year == 2021) %>%
  mutate(iso2 = tolower(ISO2),
         info = paste(country))  # Assuming you also meant to use 'inflation' here

tmap_mode("view")

tm <- tm_shape(map_for_V_2021) +
  tm_polygons("bmi_usdprice", palette = "YlOrBr", title = "BMI (US$)", id = "info", legend.is.portable = FALSE) +
  tm_facets(by = "continent", ncol = 2) + 
  tm_layout(legend.show=FALSE)+
  tm_basemap("CartoDB.DarkMatterNoLabels")

tm

2.2 Interactivity for Binned BMI

Code

library(leaflet)
library(RColorBrewer)
map_for_V_2021 <- map_for_V %>%
  filter(year == 2021)
map_for_V_2021$iso2 <- tolower(map_for_V_2021$ISO2)
map_for_V_2021 <- st_transform(map_for_V_2021, 4326)

mybins <- c(1,2,3,4,5,6,7,8)
pal <- colorBin(palette="YlOrBr", domain = map_for_V_2021$bmi_usdprice, na.color="transparent", bins=mybins)
map_for_V_2021$image_url <- paste0("https://flagcdn.com/16x12/", map_for_V_2021$iso2, ".png")

leaflet_map <- leaflet(data = map_for_V_2021) %>% 
  #addTiles() %>% 
  setView(lat=3, lng=10 , zoom=1) %>%
  #addProviderTiles(providers$CartoDB.DarkMatterNoLabels) %>%
  addProviderTiles(providers$NASAGIBS.ViirsEarthAtNight2012) %>%
  addPolygons(fillOpacity = 1, weight = 1.2, color = "#1D201F",
              fillColor = ~pal(bmi_usdprice),
              popup = ~paste('<img srcset="', image_url, '" style="width:20px;height:auto;">',
                country, "<br>",
                             "BMI USD Price: $", format(bmi_usdprice, nsmall = 2), "<br>", 
                             "BMI Local Price: $", format(bmi_localprice, nsmall = 2), "<br>", 
                             "GDP per Capita: $", format(gdp_per_capita, nsmall = 2), "<br>",
                             "Inflation(%): ", format(inflation, nsmall = 2)),
              label = ~as.character(country),  # This will be the tooltip text
              labelOptions = labelOptions(
                style = list("font-weight" = "bold", padding = "9px 12px"),
                direction = "auto",  # Automatically position the tooltip
                noHide = FALSE,  # The tooltip will hide when not hovered
                opacity = 1,  # Tooltip background opacity
                textsize = "10px"
                      ) 
)%>%
addLegend(pal=pal, values=~bmi_usdprice, opacity=0.8, title = "Big Mac Index", position = "bottomleft")


leaflet_map <- leaflet_map %>% 
  onRender("
    function(el, x) {
        var lastClickedLayer;

        this.eachLayer(function(layer) {
          layer.on('click', function(e) {
            if (lastClickedLayer) {
              lastClickedLayer.setStyle({fillColor: '#FFEDA0'});
            }
            // Change clicked layer's color to green
            layer.setStyle({fillColor: '#018749'});
            lastClickedLayer = layer;
          });
        });
      }
    ")

leaflet_map

Code

worldly<- st_read("C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/geospatial/TM_WORLD_BORDERS-0.3.shp")

Reading layer `TM_WORLD_BORDERS-0.3' from data source 
  `C:\FirGhaz\ISSS608-VAA\Take-home_Exercises\Take-Home_Ex04\data\geospatial\TM_WORLD_BORDERS-0.3.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 246 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6236
Geodetic CRS:  WGS 84

Code

worldly_projected <- st_transform(worldly, "+proj=robin")  
data <- read_csv("C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/bmi_data_2021_carto.csv")

2.3 Appreciating Trade Nodes and Edges through Cartogram

Code

map_pre <- worldly_projected%>%
  left_join(data)

clean <- data%>%
  select(iso_a3 ='ISO3',bmi_usd_price=bmi_usd_price, total_trade=total_trade, import_usd=import_usd, export_usd,)

# Join data to world map
map <- map_pre%>%
  left_join(clean)%>%
  select(-c(country, id, year, continent, iso_a3))%>%
  mutate(across(c(bmi_usd_price:population), ~ifelse(is.na(.), 0, .)))
  #drop_na(bmi_usd_price)

#st_write(map, "C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/map2021_data.csv")

ggplot(map, aes(fill=total_trade))+
  geom_sf() 

Code

dorl<-cartogram::cartogram_dorling(
  map, weight='total_trade', k = 3,
  m_weight = 1, itermax = 1000
)
# Set colors
col_world <- "#9CB4BF"
col_back <- "#1D201F"

# Set theme
theme_custom <- theme_void()+
  theme(plot.background = element_rect(fill=col_back,color=NA))

dorl<-dorl%>%
  mutate(
    # Compute area
    ar=as.numeric(st_area(dorl)),
    # Compute radius based on area
    rad=as.numeric(sqrt(ar/pi))
  )

# Extract centroids for each circle
centr <- dorl%>%
  st_centroid()%>%
  st_coordinates()

# Combine data
dorl2 <- tibble(dorl,X=centr[,1],Y=centr[,2])%>%
  arrange(-total_trade)

ggplot(map)+
  # World basemap
  geom_sf(
    worldly,mapping=aes(geometry=geometry),
    fill=col_world,color=alpha("dimgrey",0.25)
  )+
  # Draw Dorling cartogram with geom_circle()
  ggforce::geom_circle(
    data = dorl2, aes(x0 = X, y0 = Y, r = rad),
    fill=alpha("dimgrey",0.75),color=alpha("white",0.2)
  )+
  theme_custom

Code

dorl2 <- dorl2 %>%
  mutate(
    ratio_export = export_usd /total_trade,
    ratio_import = import_usd/total_trade
  )%>%
  mutate(
    rad_export=sqrt(rad*rad*ratio_export),
    rad_import=sqrt(rad*rad*ratio_import)
  )

col_export <- "#f2e901"
col_import <- "#FF4843"

ggplot(map)+
  # World basemap
  geom_sf(
    worldly,mapping=aes(geometry=geometry),
    fill=col_world,color=alpha("dimgrey",0.25)
  )+
  # Draw Dorling cartogram with geom_circle()
  ggforce::geom_circle(
    data = dorl2, aes(x0 = X, y0 = Y, r = rad),
    fill=alpha("dimgrey",0.75),color=alpha("white",0.2)
  )+
  # Draw circle for crops (or grass)
  ggforce::geom_circle(
    data = dorl2, aes(x0 = X, y0 = Y, r = rad_export),
    fill=col_export,color=NA
  )+
  theme_custom

Code

circleFun <- function(
    center=c(0,0),   # center of the circle 
    diameter=1,      # diameter 
    npoints=100,     # number of points to draw the circle
    start=0, end=2   # start point/end point
  ){
      tt <- seq(start*pi, end*pi, length.out=npoints)
      tb <- tibble(
        x = center[1] + diameter / 2 * cos(tt), 
        y = center[2] + diameter / 2 * sin(tt)
      )
    return(tb)
}

half_export <- bind_cols(
  ISO3 = rep(dorl2$ISO3[1],100),
  circleFun(
    c(dorl2$X[1],dorl2$Y[1]),dorl2$rad_export[1]*2, start=1.5, end=2.5
  ))

half_import <- bind_cols(
  ISO3 = rep(dorl2$ISO3[1],100),
  circleFun(
    c(dorl2$X[1],dorl2$Y[1]),dorl2$rad_import[1]*2, start=0.5, end=1.5
  ))

col_world <- "#9CB4BF"
col_back <- "#1D201F"

# Set theme
theme_custom <- theme_void()+
  theme(plot.background = element_rect(fill=col_back,color=NA))



# Make loop for all countries
for (i in 2:dim(dorl2)[1]){

  # Draw for exports
  temp_export <- bind_cols(
    ISO3 = rep(dorl2$ISO3[i],100),
    circleFun(
      c(dorl2$X[i],dorl2$Y[i]),dorl2$rad_export[i]*2, start=1.5, end=2.5
    ))
  # Draw for imports
  temp_import <- bind_cols(
    ISO3 = rep(dorl2$ISO3[i],100),
    circleFun(
      c(dorl2$X[i],dorl2$Y[i]),dorl2$rad_import[i]*2, start=0.5, end=1.5
    ))
  
  half_export<-half_export%>%
    bind_rows(temp_export)
  
  half_import<-half_import%>%
    bind_rows(temp_import)
}
df <- data.frame(
  x = 1:3,
  y = c(3, 2, 1),
  group = c("A", "B", "C")
)

# Make map
p <- ggplot(map)+
  # World basemap
  geom_sf(
    worldly,mapping=aes(geometry=geometry),
    fill=col_world,color=alpha("dimgrey",0.25)
    
  )+
  # Draw Dorling cartogram with geom_circle()
  ggforce::geom_circle(
    data = dorl2, aes(x0 = X, y0 = Y, r = rad),
    fill=alpha("dimgrey",0.75),color=alpha("white",0.2)
  )+
  # Draw half circle for crop with geom_polygon
  geom_polygon(
    half_export,
    mapping=aes(x,y,group=ISO3),
    fill=col_export,color=NA
  )+ 
  # Draw half circle for grass with geom_polygon
  geom_polygon(
    half_import,
    mapping=aes(x,y,group=ISO3),
    fill=col_import,color=NA
  )+ 
  theme_custom

x_range <- range(dorl2$X, na.rm = TRUE)
y_range <- range(dorl2$Y, na.rm = TRUE)
# Calculate positions for legend items based on the range
legend_x <- x_range[1] * 1.20 # Place legend slightly to the right of the X range minimum
legend_y_total_trade <- y_range[1] * 1.15 # Start legend slightly above the Y range minimum
legend_y_exports <- legend_y_total_trade + abs(y_range[2] - y_range[1]) * 0.05 # 5% above total trade
legend_y_imports <- legend_y_exports + abs(y_range[2] - y_range[1]) * 0.05

p +  annotate("point", x = legend_x, y = legend_y_total_trade, shape = 21, size = 12, fill = "dimgrey", colour = "dimgrey" ) +
  annotate("text", x = legend_x + 1, y = legend_y_total_trade, label = "Total Trade         1 billion(US)", hjust = 0.52, vjust = 1.6, colour = "grey", size=3.7, fontface = "italic", family = "Arial" ) +
  annotate("text", x = legend_x + 1, y = legend_y_exports, label = "Exports", hjust = 1.8, vjust = 1.2, colour = "#f2e901", size=3.7,fontface = "italic", family = "Arial") +
  annotate("text", x = legend_x + 1, y = legend_y_imports, label = "Imports", hjust = 1.8, vjust = 0.8, colour = "#FF4843", size=3.7, fontface = "italic", family = "Arial") +
  annotate("text", x = legend_x + 1, y = legend_y_imports, label = "Trade 2021", hjust = 1.16, vjust = -0.5, colour = "white", size=4.5, fontface = "bold", family = "Arial") +
  theme_custom

Code

bmi_node_2021 <- read.csv('C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/bmi_node_2021_add.csv')
bmi_edges_2021 <- read.csv('C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/import_export_beef_edges5-2021.csv')

The provided image is a geo-network visualization that maps the global trade network for beef, highlighting the interconnected trade relationships of over 50 countries with a detailed focus on 28 of them. The visualization uses directional arcs and varying circle sizes to depict the flow and volume of beef exports and imports between nations. Such a graphic representation allows for a quick assessment of trade balances and intensities, with larger circles likely indicating major trade hubs or high trade volumes.

This visualization is particularly valuable in the context of economic studies, including the analysis of the Big Mac Index. The Big Mac Index is an informal benchmark used to compare the purchasing power between currencies by using the price of a Big Mac as a standard good. By overlaying the beef trade data, one can derive insights into local market prices of beef, a critical component of the Big Mac, hence affecting the index’s outcomes. For instance, beef-exporting countries might have cheaper Big Mac prices due to lower local beef costs, while importing countries could face higher prices reflecting the import costs.

Overall, the map serves as a powerful tool in illustrating how trade flows influence economic indicators like the Big Mac Index. It encapsulates the complexity of global trade in a digestible format, offering economists and decision-makers a clear visualization of how commodity movements impact local market economics and consumer prices. This, in turn, can guide more informed policy decisions and economic analyses.

Code

library(sf)
library(ggplot2)
library(ggforce) 
library(maps) 
library(geosphere)

# Assuming 'bmi_node_2021' is your dataframe and it has columns 'LON' and 'LAT'
# Make sure 'bmi_node_2021' is a data frame and 'LON' and 'LAT' are not factors
bmi_node_2021a <- data.frame(bmi_node_2021)
bmi_node_2021a$LON <- as.numeric(as.character(bmi_node_2021a$LON))
bmi_node_2021a$LAT <- as.numeric(as.character(bmi_node_2021a$LAT))

# Convert your data frame to an sf object
node_sf <- st_as_sf(bmi_node_2021a, coords = c("LON", "LAT"), crs = 4326)
node_sf_projected <- st_transform(node_sf, st_crs(worldly_projected))


bmi_edges_2021$source_lon <- as.numeric(bmi_edges_2021$source_lon)
bmi_edges_2021$source_lat <- as.numeric(bmi_edges_2021$source_lat)
bmi_edges_2021$target_lon <- as.numeric(bmi_edges_2021$target_lon)
bmi_edges_2021$target_lat <- as.numeric(bmi_edges_2021$target_lat)

col_world <- "#9CB4BF"
col_back <- "#1D201F"

# Set theme
theme_custom <- theme_void()+
  theme(plot.background = element_rect(fill=col_back,color=NA))

# Create a LINESTRING for each row in bmi_edges_2021
edges_linestrings <- lapply(1:nrow(bmi_edges_2021), function(i) {
  row <- bmi_edges_2021[i, ]
  st_linestring(matrix(c(row$source_lon, row$source_lat, row$target_lon, row$target_lat), ncol = 2, byrow = TRUE))
})

# Convert the list of LINESTRINGs to an sf object
edges_sf <- st_sfc(edges_linestrings, crs = 4326)

# Convert to a data frame structure for sf, if needed
edges_sf_df <- st_sf(geometry = edges_sf)

# Transform to the same CRS as the base map
edges_sf_projected <- st_transform(edges_sf_df, st_crs(worldly_projected))



col.1 <- adjustcolor("orangered", alpha=0.04)
col.2 <- adjustcolor("yellow", alpha=0.04)
edge.pal <- colorRampPalette(c(col.1, col.2), alpha = TRUE)

# Assuming 'Weight' is normalized between 0 and 1 for color mapping
max_weight <- max(bmi_edges_2021$Value, na.rm = TRUE)
edge.col <- edge.pal(1000000)[round(1000000 * bmi_edges_2021$Value / max_weight)]

# Generate arc data
arc_data <- lapply(1:nrow(bmi_edges_2021), function(i) {
  row <- bmi_edges_2021[i,]
  arc <- gcIntermediate(
    c(row$source_lon, row$source_lat),
    c(row$target_lon, row$target_lat),
    n=300, addStartEnd=TRUE, breakAtDateLine=F
  )
  data.frame(lon = arc[,1], lat = arc[,2], weight = row$Value, edge_col = edge.col[i])
}) %>%
  do.call(rbind, .)

# Convert arc data to an sf object and project
arc_sf <- st_as_sf(arc_data, coords = c("lon", "lat"), crs = 4326)
arc_sf_projected <- st_transform(arc_sf, st_crs(worldly_projected))


q <- ggplot(map) +
  geom_sf(data = worldly_projected, fill = col_world, color = alpha("dimgrey", 0.25)) +
  geom_sf(data = node_sf_projected, aes(geometry = geometry), color = "red", size = 1) + 
  geom_sf(data = arc_sf_projected, aes(color = edge_col), size = 0.0001) +  # Assuming 'edge_col' is part of 'arc_sf_projected'
  scale_color_identity() +

  # Draw Dorling cartogram with geom_circle()
  ggforce::geom_circle(
    data = dorl2, aes(x0 = X, y0 = Y, r = rad),
    fill=alpha("dimgrey",0.75),color=alpha("white",0.2)
  )+
  # Draw half circle for crop with geom_polygon
  geom_polygon(
    half_export,
    mapping=aes(x,y,group=ISO3),
    fill=col_export,color=NA
  )+ 
  # Draw half circle for grass with geom_polygon
  geom_polygon(
    half_import,
    mapping=aes(x,y,group=ISO3),
    fill=col_import,color=NA
  )+ 
  theme_custom

x_range <- range(dorl2$X, na.rm = TRUE)
y_range <- range(dorl2$Y, na.rm = TRUE)
# Calculate positions for legend items based on the range
legend_x <- x_range[1] * 1.20 # Place legend slightly to the right of the X range minimum
legend_y_total_trade <- y_range[1] * 1.15 # Start legend slightly above the Y range minimum
legend_y_exports <- legend_y_total_trade + abs(y_range[2] - y_range[1]) * 0.05 # 5% above total trade
legend_y_imports <- legend_y_exports + abs(y_range[2] - y_range[1]) * 0.05

q + annotate("point", x = legend_x, y = legend_y_total_trade, shape = 21, size = 9, fill = "dimgrey", colour = "dimgrey") +
  annotate("text", x = legend_x + 1, y = legend_y_total_trade, label = "Total Trade       1 billion(US)", hjust = 0.49, vjust = 1.3, colour = "grey", size=3.5, fontface = "italic", family = "Arial" ) + 
  annotate("text", x = legend_x + 1, y = legend_y_total_trade, label = "Beef Trade --", hjust = 0.97, vjust = 2.6, colour = "orangered", size=3.5, fontface = "bold", family = "Arial" ) +
  annotate("text", x = legend_x + 1, y = legend_y_exports, label = "Exports", hjust = 1.8, vjust = 0.9, colour = "#f2e901", size=3.5,fontface = "italic", family = "Arial") +
  annotate("text", x = legend_x + 1, y = legend_y_imports, label = "Imports", hjust = 1.8, vjust = 0.6, colour = "#FF4843", size=3.5, fontface = "italic", family = "Arial") +
  annotate("text", x = legend_x + 1, y = legend_y_imports, label = "Trade Network 2021", hjust = 0.58, vjust = -0.6, colour = "white", size=4.3, fontface = "italic", family = "Arial") +
  
theme_custom

Code

#ggsave("geotrade1.png", plot = r, width = 10, height = 6, dpi = 300)

To build Interactivity

Code

bmi_node_2021 <- read.csv('C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/bmi_node_2021_add.csv')
bmi_edges_2021 <- read.csv('C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/data/import_export_beef_edges5-2021.csv')

img <- png::readPNG("C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/images/geotrade1.png")

img_base64 <- base64enc::dataURI(file = "C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/images/geotrade1.png", mime = "image/png")


fig <- plot_ly() %>% 
  layout(
    images = list(
      list( 
        source = img_base64,
        xref = "paper",
        yref = "paper",
        x = -0.05,
        y = 1.05, 
        sizex = 1.1,
        sizey = 1.13,#1.11
        sizing = "stretch",
        opacity = 1,
        layer = "below"
      )
    ),
    xaxis = list(range = c(-180, 180), showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
    yaxis = list(range = c(-90, 90), showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
 geo = list(
      showland = FALSE, 
      showcountries = FALSE,
      showocean = FALSE,
      showcoastlines = FALSE,
      bgcolor = 'rgba(0,0,0,0)',
      projection = list(type = 'robinson',
      lataxis = list(range = c(-90, 90), fixedrange = TRUE), 
      lonaxis = list(range = c(-180, 180), fixedrange = TRUE),
      scope = 'world',
      showframe = FALSE
)), dragmode = FALSE)

#Initialize vectors to store all edge coordinates with NA as separators
#all_lons <- c()
#all_lats <- c()
#all_texts <- c()  # For hover text

#Populate the vectors
#for(i in 1:nrow(bmi_edges_2021)) {
 # all_lons <- c(all_lons, bmi_edges_2021$source_lon[i], bmi_edges_2021$target_lon[i], NA)
  #all_lats <- c(all_lats, bmi_edges_2021$source_lat[i], bmi_edges_2021$target_lat[i], NA)
 # all_texts <- c(all_texts, paste("Origin:", bmi_edges_2021$Origin[i], 
                             #       "| Destination:", bmi_edges_2021$Partners[i]), "", NA)
#}

#edge_hover_text <- with(bmi_edges_2021, paste("", Origin, "|", Partners))
 
#fig <- fig %>% add_trace(
   #type = 'scattergeo',
 # mode = 'lines',
  # lon = all_lons,
 # lat = all_lats,
 #  text = all_texts,  # Assigning hover text
  #  hoverinfo = 'text',
  #  line = list(color = 'orange', width = 0.006),
  #  showlegend = FALSE,
   # name = "Trade Links"  # Naming the edges trace
#)

bmi_node_2021$hover_text <- paste(
  "",bmi_node_2021$country,
  "<br> BMI (USD): $", bmi_node_2021$bmi_usd_price,
  "<br> Total Trade: $", bmi_node_2021$total_trade,
  "<br> Imports: $", bmi_node_2021$import_usd,
  "<br> Exports: $", bmi_node_2021$export_usd)
 
fig <- fig %>% 
  add_trace(
    data = bmi_node_2021,
    type = 'scattergeo',
    mode = 'markers',
    text = ~hover_text,
    hoverinfo = 'text',
    hoverlabel = list(bgcolor = "white", font = list(color = "black")),
    marker = list(
      size = 0.0001, # Adjust size as needed
      color = 'grey90',
       line = list(width = 1, color = 'red'),
      projection = list(type = 'robinson')
    ),
    lon = ~LON,
    lat = ~LAT,
    dragmode = FALSE,
    rotation = list(LAT = 0, LON = 0), 
    lataxis = list(range = c(-90, 90), fixedrange = TRUE), 
    lonaxis = list(range = c(-180, 180), fixedrange = TRUE),
    scope = 'world',
    showframe = FALSE,
    dragmode = FALSE
  ) 

fig <- fig%>% config(displayModeBar = FALSE)

fig

3 Visualising the Network: Exploratory Data Analysis

Emphasis here is to give them the option to play with the nodes, layout and geom_edge functions.

Code

library(tidygraph)
library(igraph)
library(ggplot2)
library(ggraph)
library(visNetwork) 

Code

bmi2021_graph <- tbl_graph(nodes = bmi_node_2021,
                           edges = bmi_edges_2021, 
                            directed = TRUE)

Basic kk Layout + Link

Code

#this are the node attributes
V(bmi2021_graph)$bmi_usd_price <- bmi_node_2021$bmi_usd_price
V(bmi2021_graph)$gdp_per_capita<- bmi_node_2021$gdp_per_capita
V(bmi2021_graph)$inflation <- bmi_node_2021$inflation

#Layout options

ggraph(bmi2021_graph, layout = "kk") + #layout options = kk, fr, nicely,
  geom_node_point(aes(fill = V(bmi2021_graph)$bmi_usd_price , size = V(bmi2021_graph)$bmi_usd_price, color = continent, alpha = 0.7), show.legend=FALSE)+
  geom_edge_link2(aes(color = Trade, width = Value), alpha = 0.4, show.legend=FALSE) + 
  scale_edge_color_manual(values = c("Import Quantity" = "#FF4843", "Export Quantity" = "#f2e901")) +
  scale_edge_width(range = c(0.5, 7)) +
  geom_node_text(aes(label = country), repel = TRUE, size = 3, hjust = 0, vjust = 0, nudge_y=-0.2, nudge_x=-0.2, color = "white") +
  theme_graph(background = '#1D201F', text_colour = 'white')

Linear Arc Layout + arc

Code

#this are the node attributes
V(bmi2021_graph)$bmi_usd_price <- bmi_node_2021$bmi_usd_price
V(bmi2021_graph)$gdp_per_capita<- bmi_node_2021$gdp_per_capita
V(bmi2021_graph)$inflation <- bmi_node_2021$inflation

#Layout options

ggraph(bmi2021_graph, layout = "linear") + #layout options = kk, fr, nicely,
  geom_node_point(aes(fill = V(bmi2021_graph)$bmi_usd_price , size = V(bmi2021_graph)$bmi_usd_price, colour = continent, alpha = 0.7), show.legend=FALSE)+
  geom_edge_arc(aes(color = Trade, width = Value), alpha = 0.4, show.legend=FALSE) + 
  scale_edge_color_manual(values = c("Import Quantity" = "#FF4843", "Export Quantity" = "#f2e901")) +
  scale_edge_width(range = c(0.5, 7)) +
  geom_node_text(aes(label = country), repel = FALSE, size = 3.4, hjust = 0, vjust = 0, angle = 270, color = "white", nudge_y=-3, nudge_x=-0.3) +
  theme_graph(background = '#1D201F', text_colour = 'white')+
  theme(legend.position="none",
    plot.margin=unit(c(0,0,0.0,0), "null"),
    panel.spacing=unit(c(0,0,3.4,0), "null")
  ) 

Circular Arc linear + arc

Code

node_positions <- as.data.frame(create_layout(bmi2021_graph, layout = "linear", circular = TRUE))

# Calculate angles for each node position
node_positions$angle <- atan2(node_positions$y, node_positions$x) * (180 / pi)

# Adjust angle for text alignment
node_positions$label_angle <- ifelse(node_positions$angle > 90 | node_positions$angle < -90, node_positions$angle + 180, node_positions$angle)

# Define if the text should be above or below the node based on its angle
node_positions$hjust <- ifelse(node_positions$angle > 90 | node_positions$angle < -90, 1, 0)

#this are the node attributes
V(bmi2021_graph)$bmi_usd_price <- bmi_node_2021$bmi_usd_price
V(bmi2021_graph)$gdp_per_capita<- bmi_node_2021$gdp_per_capita
V(bmi2021_graph)$inflation <- bmi_node_2021$inflation
#Layout options

p <- ggraph(bmi2021_graph, layout = "linear", circular = TRUE) + #layout options = kk, fr, nicely,
  geom_node_point(aes(fill = V(bmi2021_graph)$bmi_usd_price , size = V(bmi2021_graph)$bmi_usd_price, colour = continent, alpha = 0.7), show.legend=FALSE)+
  geom_edge_arc(aes(color = Trade, width = Value), alpha = 0.3, show.legend=FALSE) + 
  scale_edge_color_manual(values = c("Import Quantity" = "#FF4843", "Export Quantity" = "#f2e901")) +
  scale_edge_width(range = c(0.5, 8)) +
   geom_node_text(data = node_positions, aes(x = x*1.04, y = y*1.04, label = country, angle = label_angle, hjust = hjust), color = "white", size = 2.65) +
  theme_graph(background = '#1D201F', text_colour = 'white')+
  theme(legend.position="none",
    plot.margin=unit(c(0.0,0.0,0,0), "null"),
    panel.spacing=unit(c(0,0,0,0,0), "null")
  ) +
    expand_limits(x = c(-2.3,2.3), y = c(-1.7, 1.4))
  #guides(color = guide_legend(title = "Continents") +
  #edge_width = guide_legend(override.aes = list(color = "white"))) 
 #coord_fixed() 
p

Layout Grid & Fan

Code

#this are the node attributes
V(bmi2021_graph)$bmi_usd_price <- bmi_node_2021$bmi_usd_price
V(bmi2021_graph)$gdp_per_capita<- bmi_node_2021$gdp_per_capita
V(bmi2021_graph)$inflation <- bmi_node_2021$inflation

#Layout options

ggraph(bmi2021_graph, layout = "grid") + #layout options = kk, fr, nicely,
  geom_node_point(aes(fill = V(bmi2021_graph)$bmi_usd_price , size = V(bmi2021_graph)$bmi_usd_price, color = continent, alpha = 0.7), show.legend=FALSE)+
  geom_edge_fan(aes(color = Trade, width = Value), alpha = 0.4, show.legend=FALSE) + 
  scale_edge_color_manual(values = c("Import Quantity" = "#FF4843", "Export Quantity" = "#f2e901")) +
  scale_edge_width(range = c(0.5, 5)) +
  geom_node_text(aes(label = country), repel = TRUE, size = 3, hjust = 0, vjust = 0, nudge_y=-0.02, nudge_x=0, color = "white") +
  theme_graph(background = '#1D201F', text_colour = 'white')

Code

  #guides(color = guide_legend(title = "Continents") +
  #edge_width = guide_legend(override.aes = list(color = "white"))) 
 #coord_fixed() 
#

Code

library(visNetwork)
library(htmlwidgets)

nodes_data <- data.frame(id = V(bmi2021_graph)$country, 
                         label = V(bmi2021_graph)$country,
                         value = V(bmi2021_graph)$bmi_usd_price, #to change to inflation, gdp
                         group = V(bmi2021_graph)$continent,
                         title = paste("<b>BMI:</b>", V(bmi2021_graph)$bmi_usd_price,
                                       "<br><b>GDP:</b>", V(bmi2021_graph)$gdp_per_capita,
                                       "<br><b>Inf:</b>", V(bmi2021_graph)$inflation))

edges_data <- data.frame(from = bmi_edges_2021$Origin, 
                         to = bmi_edges_2021$Partners,
                         value = bmi_edges_2021$Value,
                         title = bmi_edges_2021$Trade,
                         arrows = 'to')

network <- visNetwork(nodes_data, edges_data, width = "100%", height = "400px") %>%
  visNodes(size = V(bmi2021_graph)$bmi_usd_price, shadow = TRUE,
           color = list(background = "white", border = "#2B7CE9"),
           font = list(color = "white")) %>%
  visEdges(smooth = FALSE, color = list(color = "#A5ABB6", highlight = "#FF4F09")) %>%
  visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE) %>%
  visLegend() %>%
  visInteraction(navigationButtons = TRUE) %>%
  visLayout(randomSeed = 123)


network <- network %>% 
  htmlwidgets::onRender("
    function(el, x) {
      el.style.backgroundColor = 'black';
    }
  ")

# Print the network to display it
network

Just for the Fun of it!

Code

library(png)
library(visNetwork)
library(htmlwidgets)
library(base64enc)

# Load the image file
img_path <- "C:/FirGhaz/ISSS608-VAA/Take-home_Exercises/Take-Home_Ex04/images/Bigmac-removebg-preview.png" 
img <- readPNG(img_path)

# Convert the image to base64
img_base64 <- base64enc::dataURI(file = img_path, mime = "image/png")

# Update the 'nodes_data' dataframe to include the image for each node
nodes_data$id <- V(bmi2021_graph)$country
nodes_data$label <- V(bmi2021_graph)$country
nodes_data$shape <- "image"
nodes_data$image <- img_base64 # Assign the same image to all nodes for simplicity

# Build the visNetwork object as usual
network <- visNetwork(nodes_data, edges_data, width = "100%", height = "400px") %>%
  visNodes(size = V(bmi2021_graph)$bmi_usd_price, shape = "image", image = nodes_data$image, borderWidth = 1, color = list(background = "white", border = "#2B7CE9"),
           font = list(color = "white")) %>%
  visEdges(smooth = FALSE, color = list(color = "#f2e901", highlight = "#FF4F09")) %>%
  visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE) %>%
  visInteraction(navigationButtons = TRUE) %>%
  visLayout(randomSeed = 123)

# Customize the background color with custom CSS
network <- network %>% 
  htmlwidgets::onRender("
    function(el, x) {
      el.style.backgroundColor = 'black';
    }
  ")

network

3.1 Network Centrality and Community Clustering

For our scope of the project, we will select 50+ countries with (1) significant trade activities and (2) complete data together with BMI(Index). A Trade Connectivity Index (TCI) will be calculated for each edge in the network will be calculated as follows:

Edge Weight Formula: Trade Connectivity Index (TCI)

This formula reflects the proportion of beef trade (exports and imports) in relation to the total trade volume between two countries. It highlights the significance of beef trade in the bilateral trade relationship, making it particularly relevant for analysing the trade networks with respect to the Big Mac Index. Hence, this will be applied to our network edges.

• Beef Exports Value and Beef Imports Value refer to the value of beef traded between the two countries forming an edge in the network.

• Total Exports Value and Total Imports Value encompass the total trade volume between these countries, providing a base to understand the relative importance of beef trade.

Node Attributes

For the nodes, you can maintain attributes that reflect each country’s economic and demographic scale, as well as its relevance to the Big Mac Index:

• BMI Index: Directly use the Big Mac Index as a node attribute to represent the pricing level of a Big Mac, serving as a proxy for purchasing power parity.

• Trade Value: This could be represented by either Total Exports or Total Imports to reflect the country’s engagement in global trade. Alternatively, you might use a composite measure such as GDP to reflect overall economic size and capacity.

We will compute the centrality and community clusters to determine their community clusters

In our exploration of global trade dynamics, specifically through the lens of beef exports and imports, we aim to uncover insights into the Big Mac Index, a whimsical gauge of purchasing power parity devised by The Economist. Analyzing beef trade patterns grants us a direct line of inquiry into the pricing mechanisms of the Big Mac, considering beef’s pivotal role in its composition. This focus allows us to parse the supply chain intricacies and market conditions that dictate the variability in Big Mac prices across different regions. The rationale behind this targeted analysis is to use beef trade as a tangible metric for understanding broader economic trends and pricing pressures that ultimately influence consumer goods’ costs.

To navigate this complex interplay, we’ve adopted a formula integrating critical trade data and economic indicators to model their relationship with the Big Mac Index. This formula captures the essence of trade balance, economic scale (granular and macro at the same time), and the specific impact of beef trade. It correlates net exports relative to GDP per capita with the proportions of beef exports and imports within total trade, providing a nuanced perspective on how trade dynamics influence market conditions and, by extension, Big Mac prices.

3.2 Computing Centrality

Code

library(ggraph)
library(igraph)


bmi2021_graph <- igraph::set_vertex_attr(bmi2021_graph, "betweenness", value = igraph::betweenness(bmi2021_graph))


g <- ggraph(bmi2021_graph, layout = "fr") +
  geom_edge_link(aes(width = Value), alpha = 0.5, show.legend = FALSE) +
  scale_edge_width(range = c(0.1, 5)) +
  geom_node_point(aes(size = betweenness), color = "black") +
  geom_node_text(aes(label = country), repel = TRUE, size = 2, hjust = 0.6, vjust = -1.8, color = "black") +
  theme_graph()


g <- g + geom_node_point(aes(color = continent))


g

Degree Centrality Degree centrality measures the number of edges connected to a node. In directed networks, you can distinguish between in-degree and out-degree.

Degree centrality is a measure used in network analysis to quantify the importance or influence of a particular node within a network. It is based on the number of connections, or edges, that a node has to other nodes. The central concept behind degree centrality is simple: nodes with more connections are considered more central and potentially more influential within the network.

There are two types of degree centrality:

In-Degree Centrality: This measures the number of incoming connections to a node. It can be particularly relevant in directed networks where the direction of the connection matters. A high in-degree centrality indicates that a node is a major target within the network, receiving many connections from other nodes. This can signify a node of high interest or popularity.

Out-Degree Centrality: This measures the number of outgoing connections from a node. Like in-degree centrality, it is applicable in directed networks. A high out-degree centrality signifies that a node actively reaches out to many other nodes, which can indicate a source or distributor of information, goods, or influence within the network

Code

bmi2021_graph <- igraph::set_vertex_attr(bmi2021_graph, "degree", value = igraph::degree(bmi2021_graph))

# Plotting the graph with ggraph
g <- ggraph(bmi2021_graph, layout = "fr") +
  geom_edge_link(aes(width = Value), alpha = 0.5,show.legend = FALSE) + 
  scale_edge_width(range = c(0.1, 5)) +
  geom_node_point(aes(size = degree, color = continent)) + # Size nodes by degree centrality and color by continent
  geom_node_text(aes(label = country), repel = TRUE, size = 2, hjust = 0.5, color = "black", fontface = "bold") +
  theme_graph() +
  scale_size_continuous(range = c(1, 10)) + # Adjust the range for the size of the nodes
  labs(size = "Degree Centrality") # Label for the size scale

# Display the graph
g

Degree centrality is a straightforward but powerful concept in network analysis, useful for identifying key nodes that might play critical roles in the dissemination of information, disease transmission, social network influence, and more within a network.

Eigenvector Centrality Eigenvector centrality measures a node’s influence based on the principle that connections to high-scoring nodes contribute more to the score of the node in question.

Code

eigenvector_centrality <- eigen_centrality(bmi2021_graph)$vector

print(summary(eigenvector_centrality))

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0000108 0.0010364 0.0061979 0.0611162 0.0231699 1.0000000 

Code

V(bmi2021_graph)$eigenvector <- eigenvector_centrality

print(head(V(bmi2021_graph)$eigenvector))

[1] 0.1068683 1.0000000 0.6944934 0.3322083 0.2881434 0.1124052

Code

# Now, use ggraph to visualize the network
g_eigen <- ggraph(bmi2021_graph, layout = "fr") + 
  geom_edge_link(aes(width = Value), alpha = 0.5,show.legend = FALSE) +
  scale_edge_width(range = c(0.1, 5)) +
  geom_node_point(aes(size = eigenvector, color = continent)) +
  geom_node_text(aes(label = country), repel = TRUE, size = 3, color = "black") +
  theme_graph() +
  scale_size_continuous(name = "Eigenvector Centrality") +
  guides(color = guide_legend("Continent"))

# Print the graph with eigenvector centrality
print(g_eigen)

Eigenvector centrality is useful for identifying influential nodes in a network where connections to high-scoring nodes contribute more to the score of the node than equal connections to low-scoring nodes.

Closeness centrality

Putting the results in a dataframe.

                   node degree betweenness   closeness  eigenvector
1                Sweden     56  51.0000000 0.004424779 1.000000e+00
2                Israel     66 100.0000000 0.005714286 6.944934e-01
3                Canada     52  24.0000000 0.250000000 3.322083e-01
4               Germany     44  90.8183669 0.007812500 2.881434e-01
5              Malaysia     22   1.5000000 0.142857143 1.500285e-01
6           Netherlands     19  38.1816331 0.012500000 1.124052e-01
7           Switzerland     32   1.5000000 0.003597122 1.068683e-01
8               Belgium     62 173.4975715 0.023255814 1.005888e-01
9                France     29  42.2459815 0.009090909 8.165931e-02
10        United States     51  71.7345151 0.052631579 4.604534e-02
11            Hong Kong      6   1.5000000 0.142857143 4.091687e-02
12              Austria     41  63.1471309 0.012820513 3.472823e-02
13            Indonesia      4   1.5000000 0.142857143 2.727791e-02
14       United Kingdom     40  95.8733867 0.017857143 2.316993e-02
15            Singapore     14   8.1621552 0.014285714 2.145370e-02
16                Spain     18  25.6956297 0.014492754 2.078596e-02
17              Denmark     16  31.1744346 0.090909091 1.970763e-02
18            Argentina     20   6.9444526 0.066666667 1.797603e-02
19                India      2   1.5000000 0.142857143 1.363896e-02
20            Australia     16  23.1709300 0.090909091 1.071837e-02
21             Portugal     16   9.0733029 0.033333333 1.053569e-02
22              Lebanon     35  70.3148125 0.021739130 9.747172e-03
23               Greece     16   9.1393395 0.125000000 9.049026e-03
24           Costa Rica     61 133.3493109 0.200000000 7.821428e-03
25                Japan     11  27.7088178 0.066666667 7.008999e-03
26                Qatar     14   1.5000000 0.013333333 6.658762e-03
27              Uruguay     25  21.3652549 0.025000000 6.197949e-03
28          New Zealand     10  32.4385155 0.043478261 6.096445e-03
29               Kuwait      7  19.9409967 0.333333333 5.675541e-03
30             Thailand     23  23.2725229 0.500000000 3.172497e-03
31       Czech Republic     34  26.4117937 0.142857143 2.833332e-03
32               Turkey     16   1.5000000 0.111111111 2.568897e-03
33               Poland      4   1.5000000 0.013333333 1.902503e-03
34                China      4   1.5000000 0.013333333 1.902503e-03
35               Brazil     17  13.3362885 0.250000000 1.837269e-03
36          South Korea     10  11.2371524 0.038461538 1.823930e-03
37         Saudi Arabia     14  16.0000000 0.033333333 1.774966e-03
38              Ukraine      8   1.5000000 0.111111111 1.284449e-03
39                Ghana      7   1.5000000 0.111111111 1.123893e-03
40              Hungary     20  16.0000000 0.100000000 1.036376e-03
41 United Arab Emirates      5   1.5000000 0.000100000 1.000436e-03
42         South Africa      6   1.5000000 0.111111111 9.633364e-04
43                 Peru      8   2.0000000 0.333333333 8.962643e-04
44               Jordan     14  14.0000000 0.333333333 8.769263e-04
45               Russia      4  30.0025813 0.111111111 8.339978e-04
46              Croatia      6   1.5000000 0.020408163 7.633785e-04
47                Chile     14  26.0404078 0.055555556 2.480581e-04
48                Egypt      6   1.5000000 0.166666667 2.262896e-04
49              Bahrain      5   0.7227148 1.000000000 2.173124e-04
50               Mexico      7  12.0000000 1.000000000 2.097795e-04
51                 Oman     10   1.5000000 0.000100000 3.243397e-05
52          Philippines      4   1.5000000 1.000000000 1.594200e-05
53             Viet Nam      3   1.5000000 1.000000000 1.084993e-05

Degree Centrality: The simplest form of centrality, degree centrality, counts how many connections (edges) a node has. In directed networks, you can further distinguish between in-degree (number of incoming edges) and out-degree (number of outgoing edges). Nodes with a high degree centrality are considered highly connected or active within the network.

Closeness Centrality: This measure calculates how close a node is to all other nodes in the network by considering the shortest paths. A node with a high closeness centrality can quickly interact with all others because it has the shortest average distance to all other nodes.

Betweenness Centrality: Betweenness centrality quantifies the number of times a node acts as a bridge along the shortest path between two other nodes. Nodes with high betweenness centrality have significant control over information flow within the network because they lie on many shortest paths between other nodes.

Eigenvector Centrality: This measure considers not just the quantity of connections a node has, but also the quality. A node scores high in eigenvector centrality if it is connected to other nodes that are themselves central within the network. This measure reflects the intuition that not all connections are equal, and being connected to highly connected nodes contributes more to a node’s score.

Code

library(dplyr)

# Assuming `centrality_df` is your dataframe with centrality measures
library(dplyr)
library(tidyr)


# Assuming `centrality_measures` is your dataframe with centrality measures
normalized_centrality_measures <- centrality_measures %>%
  # Replace NaN values in closeness with 0 (or another appropriate value)
  mutate(closeness = ifelse(is.nan(closeness), 0.0001, closeness)) %>%
  mutate(betweenness = case_when(betweenness == 0 ~ 1.5, TRUE ~ betweenness)) %>%
  #mutate(eigenvector = case_when(eigenvector == 0 ~ 1.5, TRUE ~ eigenvector))
  mutate(across(c(degree, betweenness, closeness, eigenvector),
                ~(. - min(.)) / (max(.) - min(.))))

Code

library(ggplot2)
library(tidyr)
library(ggplot2)
library(tidyr)
library(dplyr)

# Assuming `normalized_centrality_measures` is your original data frame
selected_countries <- c("Argentina", "Australia", "Brazil", "United Kingdom", 
                        "Canada", "Chile", "China", "Czech Republic", 
                        "Denmark", "Hong Kong", "Hungary", "Indonesia", 
                        "Japan", "Malaysia", "Mexico", "New Zealand", 
                        "Peru", "Philippines", "Poland", "Russia", 
                        "Singapore", "South Africa", "South Korea", 
                        "Sweden", "Switzerland", "Thailand", "Turkey", 
                        "United States")

# Pivot the data to long format
long_centrality_measures <- normalized_centrality_measures %>%
  filter(node %in% selected_countries) %>%
  pivot_longer(cols = c(degree, betweenness, closeness, eigenvector),
               names_to = "measure",
               values_to = "value")

# Create a function to generate plots for a given measure
generate_plot <- function(measure) {
  plot <- ggplot(long_centrality_measures %>% filter(measure == !!measure), 
                 aes(x = reorder(node, value), y = value)) +
    geom_col(aes(fill = value)) + 
    scale_fill_gradient(low = "#f2e901", high = "#FF2F09") +
    
    labs(x = "Node", y = paste(measure, "Centrality (Normalized)"), 
         title = paste(measure, "Centrality of Nodes")) +
    coord_flip() + 
    theme_minimal() +
    theme(
      axis.text.y = element_text(size = 6),  # Adjust the y-axis font size
      strip.text.y = element_text(size = 6),
      legend.position = "none"# Adjust the facet label font size if needed
    ) +
    scale_y_continuous(limits = c(0, 1))  
  return(plot)
}

# Generate the plots
degree_plot <- generate_plot("degree")
betweenness_plot <- generate_plot("betweenness")
closeness_plot <- generate_plot("closeness")
eigenvector_plot <- generate_plot("eigenvector")

# Print the plots
degree_plot

Code

betweenness_plot

Code

closeness_plot

Code

eigenvector_plot

Code

selected_countries <- c("Argentina", "Australia", "Brazil", "United Kingdom", 
                        "Canada", "Chile", "China", "Czech Republic", 
                        "Denmark", "Hong Kong", "Hungary", "Indonesia", 
                        "Japan", "Malaysia", "Mexico", "New Zealand", 
                        "Peru", "Philippines", "Poland", "Russia", 
                        "Singapore", "South Africa", "South Korea", 
                        "Sweden", "Switzerland", "Thailand", "Turkey", 
                        "United States")

long_centrality_measures <- normalized_centrality_measures %>%
  filter(node %in% selected_countries) %>%
  pivot_longer(cols = c(degree, betweenness, closeness, eigenvector),
               names_to = "measure",
               values_to = "value")

###################################################################################

# Merge back the original degree centrality measure for hover info
long_centrality_measures_with_original <- long_centrality_measures %>%
  left_join(centrality_measures %>% select(node, original_degree = degree), by = "node")

# Generate the plot for degree measure
degree_data <- long_centrality_measures_with_original %>%
  filter(measure == "degree")

degree_plot <- ggplot(degree_data, aes(x = reorder(node, value), y = value, text = paste("Country:", node, "<br>Degree Centrality:", original_degree))) +
  geom_col(aes(fill = value)) +
  scale_fill_gradient(low = "#f2e901", high = "#FF2F09") +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Node", y = "Degree Centrality (Normalized)", title = "Degree Centrality of Nodes") +
  coord_flip() +
  theme_minimal() +
  theme(
      axis.text.y = element_text(size = 5),  # Adjust the y-axis font size
      strip.text.y = element_text(size = 5),
      legend.position = "none"# Adjust the facet label font size if needed
    )

# Convert to Plotly for interactivity
degree_plotly <- ggplotly(degree_plot, tooltip = "text")

# Print the interactive plot
degree_plotly

Code

# Merge back the original closeness centrality measure for hover info
long_centrality_measures_with_original <- long_centrality_measures %>%
  left_join(centrality_measures %>% select(node, original_closeness = closeness), by = "node")

# Generate the plot for closeness measure
closeness_data <- long_centrality_measures_with_original %>%
  filter(measure == "closeness")

closeness_plot <- ggplot(closeness_data, aes(x = reorder(node, value), y = value, text = paste("Country:", node, "<br> Closeness Centrality:", original_closeness))) +
  geom_col(aes(fill = value)) +
  scale_fill_gradient(low = "#f2e901", high = "#FF2F09") +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Node", y = "Closeness Centrality (Normalized)", title = "Closeness Centrality of Nodes") +
  coord_flip() +
  theme_minimal() +
  theme(
      axis.text.y = element_text(size = 5),  # Adjust the y-axis font size
      strip.text.y = element_text(size = 5),
      legend.position = "none"# Adjust the facet label font size if needed
    )

# Convert to Plotly for interactivity
closeness_plotly <- ggplotly(closeness_plot, tooltip = "text")

# Print the interactive plot
closeness_plotly

Code

long_centrality_measures_with_original <- long_centrality_measures %>%
  left_join(centrality_measures %>% select(node, original_betweenness = betweenness), by = "node")

# Generate the plot for closeness measure
betweenness_data <- long_centrality_measures_with_original %>%
  filter(measure == "betweenness")

betweenness_plot <- ggplot(betweenness_data, aes(x = reorder(node, value), y = value, text = paste("Country:", node, "<br> Betweenness Centrality:", original_betweenness))) +
  geom_col(aes(fill = value)) +
  scale_fill_gradient(low = "#f2e901", high = "#FF2F09") +
  scale_y_continuous(limits = c(0, 1)) +
  labs(x = "Node", y = "Betweenness Centrality (Normalized)", title = "betweenness Centrality of Nodes") +
  coord_flip() +
  theme_minimal() + theme(
      axis.text.y = element_text(size = 5),  # Adjust the y-axis font size
      strip.text.y = element_text(size = 5),
      legend.position = "none"# Adjust the facet label font size if needed
    )

# Convert to Plotly for interactivity
betweenness_plotly <- ggplotly(betweenness_plot, tooltip = "text")

# Print the interactive plot
betweenness_plotly

Code

long_centrality_measures_with_original <- long_centrality_measures %>%
  left_join(centrality_measures %>% select(node, original_eigenvector = eigenvector), by = "node")

# Generate the plot for eigenvector measure
eigenvector_data <- long_centrality_measures_with_original %>%
  filter(measure == "eigenvector")

eigenvector_plot <- ggplot(eigenvector_data, aes(x = reorder(node, value), y = value, text = paste("Country:", node, "<br>Eigenvector Centrality:", original_eigenvector))) +
  geom_col(aes(fill = value)) +
  scale_fill_gradient(low = "#f2e901", high = "#FF2F09") +
  labs(x = "Node", y = "Eigenvector Centrality (Normalized)", title = "Centrality of Nodes") +
  coord_flip() +
  theme_minimal() + theme(
      axis.text.y = element_text(size = 5),  # Adjust the y-axis font size
      strip.text.y = element_text(size = 5),
      legend.position = "none"# Adjust the facet label font size if needed
    )

# Convert to Plotly for interactivity
eigenvector_plotly <- ggplotly(eigenvector_plot, tooltip = "text")

# Print the interactive plot
eigenvector_plotly

Code

combined_plot <- subplot(degree_plotly, betweenness_plotly,closeness_plotly, eigenvector_plotly, nrows = 2, margin = 0.05)

# Print the combined plot
combined_plot

In the Shiny App, we will allow the audience to explore these statistical methods to gain insights from the networks.

Betweenness	Degree	EigenVector	Closeness
Nodes like the United Kingdom, United States and Sweden have high betweenness centrality, suggesting that they act as significant connectors or bridges within the network, potentially controlling the flow of trade.	Nodes such as the Sweden, Canada, United States and UK have high degrees, meaning they have numerous trade connections, possibly making them central hubs in the beef trade network.	The eigenvector centrality considers not just the number but the quality of connections. Nodes with high eigenvector centrality, like the Sweden and Canada, are connected to other well-connected nodes, hinting at influential trade cliques.	High closeness centrality for nodes like the Philippines, Thailand and Mexico suggests that they can quickly interact or trade with all other nodes, indicating efficiency in their trade operations.
Some countries much higher centrality, which could mean the trade network relies heavily on these nodes, possibly creating vulnerability to disruptions.	Degree centrality can suggest potential market influence, with highly connected nodes being able to leverage their position for competitive advantage.	This measure reflects the potential for a node to access and influence the broader network through its connections.	Closeness centrality offers insight into the speed at which a country can react to supply and demand changes across the global network
Countries with lower betweenness centrality might have more direct trade links or may not be critical transit points within the network.	A spread of degree centrality across the network would indicate a democratized trading environment, but a concentration suggests a few key players dominate.	Low eigenvector centrality could indicate peripheral or less influential roles within the trading network, possibly representing newer or more specialized markets.	Low closeness centrality might point to potential delays in trade flows or inefficiencies, possibly due to geographical or logistical factor

:::callout: Summary of Insights 1. Network Roles and Big Mac Pricing:

Countries central in the beef trade network, indicated by high betweenness or degree centrality, may have more stable and competitive beef pricing due to their numerous trade links, which could affect their local Big Mac prices. The presence of key influencers, as suggested by eigenvector centrality, could indicate the potential for price setting or market influence, which can trickle down to the pricing of beef and related products like the Big Mac.

2. Trade Efficiency and Cost:

Closeness centrality can signal efficient trade practices and the potential for rapid adjustment to market changes, which could lead to more competitive and stable pricing for beef, affecting the Big Mac Index. Any delays or inefficiencies in beef imports could increase costs for local producers, potentially raising Big Mac prices in countries with lower closeness centrality.

3. Cluster Analysis and Market Groupings:

Clustering helps identify groups of countries with similar trade patterns, which can show regional variations in Big Mac prices due to similarities in trade dynamics and economic profiles. Understanding these clusters can provide insights into regional pricing strategies for the Big Mac and indicate how economic or trade changes in one country might impact others within the same cluster. :::

3.3 Computing Community Indices

Tidygraph package inherits many of the community detection algorithms imbedded into igraph We will utilise 2 community detection algo: (1) Walktrap and (2) Spinglass (group_spinglass).Some community algorithms are designed to take into account direction or weight, while others ignore it.

WalkTrap The Walktrap algorithm is another method for detecting communities in graphs. It attempts to find densely connected subgraphs (communities) in a graph based on random walks. The idea is that short random walks tend to stay in the same community.

Code

# Calibrating Parameters for Walktrap
# Choosing the weights
#E(bmi2021_graph)$Value
#E(bmi2021_graph)$betweenness 
#E(bmi2021_graph)$eigenvector 
#E(bmi2021_graph)$closenesss
#E(bmi2021_graph)$degree 

weights = E(bmi2021_graph)$Value
steps = 30 # (2 to 500)

walktrap_communities <- cluster_walktrap(bmi2021_graph, 
                                         steps = steps,
                                         weights = weights)
# Assign community membership to nodes for visualization
V(bmi2021_graph)$community <- as.factor(membership(walktrap_communities))

g <- ggraph(bmi2021_graph, layout = 'kk') +
  geom_edge_link2(aes(width = Value), alpha = 0.5) +
  scale_edge_width(range = c(0.5, 8)) +
  geom_node_point(aes(fill = bmi_usd_price, size = bmi_usd_price, color = community), show.legend = FALSE) +
  geom_node_text(aes(label = country), repel = TRUE, size = 2, hjust = 0.6, vjust = -1.8, color = "black") +
  theme_graph() +
  labs(fill = "BMI USD Price", color = "Community")  

g

This code performs a sweep over a range of steps values for the Walktrap community detection algorithm, calculates modularity scores for each, and visualizes these scores to identify the optimal steps value for community detection in your graph. The ggplot2 visualization will help you easily spot the steps value that maximizes the modularity score, indicating an optimal community structure.

Code

# Assuming 'bmi2021_graph' is your graph object and 'weights' are the edge weights

# Define the range for steps
steps_range <- seq(2, 500, by = 50) # Adjust the step for finer or coarser granularity

# Initialize a vector to store modularity scores
modularity_scores <- numeric(length(steps_range))

# Sweep over steps
for (i in seq_along(steps_range)) {
  steps <- steps_range[i]
  # Walktrap community detection
  walktrap_communities <- cluster_walktrap(bmi2021_graph, steps = steps, weights = E(bmi2021_graph)$Value)
  # Calculate modularity
  mod_score <- modularity(walktrap_communities, weights = E(bmi2021_graph)$Value)
  # Store the modularity score
  modularity_scores[i] <- mod_score
}

# Visualization
library(ggplot2)

# Create a data frame for plotting
modularity_df <- data.frame(Steps = steps_range, Modularity = modularity_scores)

# Line plot of modularity scores vs. steps
ggplot(modularity_df, aes(x = Steps, y = Modularity)) +
  geom_line() + 
  geom_point(shape = 1) + # Adds points to the line
  theme_minimal() +
  labs(x = "Steps", y = "Modularity", title = "Effect of Steps on Modularity in Walktrap Community Detection")

Code

# Optionally, you may want to highlight the steps value that leads to the highest modularity score
max_modularity_steps <- steps_range[which.max(modularity_scores)]
max_modularity <- max(modularity_scores)

# Print out the optimal steps value and its modularity score for reference
cat("Optimal Steps:", max_modularity_steps, "- Max Modularity:", max_modularity, "\n")

Optimal Steps: 52 - Max Modularity: 0.6343335 

Calling out the communities - These communities will be passed on for regression.

Code

community_membership <- membership(walktrap_communities)
node_names <- V(bmi2021_graph)$country 

if(length(community_membership) == length(node_names)) {
    node_stats <- data.frame(
        node = node_names,
        community = as.factor(community_membership)
    )
  
    
    countries_in_community <- aggregate(node ~ community, data = node_stats, FUN = function(x) paste(x, collapse = ", "))
  
    # Print the list of countries by community
    print(countries_in_community)
} else {
    print("Mismatch in lengths of community memberships and node names.")
}

  community
1         1
2         2
3         3
4         4
5         5
6         6
7         7
8         8
9         9
                                                                                                                                                             node
1                                                                                                                 Germany, Netherlands, Spain, Belgium, Singapore
2 Lebanon, Thailand, Czech Republic, South Korea, Chile, United Arab Emirates, Brazil, Bahrain, Peru, Jordan, Hungary, Philippines, Viet Nam, Oman, Egypt, Mexico
3                                                                                            Australia, Denmark, Costa Rica, Ukraine, South Africa, Ghana, Turkey
4                                            France, Austria, Greece, Portugal, New Zealand, Uruguay, United Kingdom, Saudi Arabia, Croatia, Qatar, Poland, China
5                                                                                                           Israel, Canada, Hong Kong, India, Malaysia, Indonesia
6                                                                                                                         United States, Kuwait, Argentina, Japan
7                                                                                                                                                     Switzerland
8                                                                                                                                                          Sweden
9                                                                                                                                                          Russia

Code

communities_walktrap <- cluster_walktrap(bmi2021_graph, weights = E(bmi2021_graph)$Value)

# Add community membership directly into the graph object
V(bmi2021_graph)$community <- as.factor(membership(communities_walktrap))

# Create nodes data frame
nodes_data <- data.frame(
  id = V(bmi2021_graph)$country,  # Assuming 'name' attribute contains the country names
  label = V(bmi2021_graph)$country,
  value = V(bmi2021_graph)$bmi_usd_price, # Adjust attribute names as necessary
  title = paste("<b>BMI:</b>", V(bmi2021_graph)$bmi_usd_price,
                "<br><b>GDP:</b>", V(bmi2021_graph)$gdp_per_capita,
                "<br><b>Inflation:</b>", V(bmi2021_graph)$inflation),
  community = V(bmi2021_graph)$community
)

# Assign colors to communities
num_communities <- length(unique(nodes_data$community))
color_palette <- colorRampPalette(RColorBrewer::brewer.pal(min(9, num_communities), "Set1"))
community_colors <- color_palette(num_communities)
color_mapping <- setNames(community_colors, levels(nodes_data$community))
nodes_data$color <- color_mapping[nodes_data$community]


edges_data <- data.frame(from = bmi_edges_2021$Origin, 
                         to = bmi_edges_2021$Partners,
                         value = bmi_edges_2021$Value,
                         title = bmi_edges_2021$Trade,
                         arrows = 'to')

# Create the network visualization
network <- visNetwork(nodes_data, edges_data) %>%
  visNodes(aes(fill = "community", color = "community"), shape = "dot", scaling = list(label = list(enabled = TRUE)), color = list(background = "white", border = "#2B7CE9"),
           font = list(color = "white")) %>%
  visEdges(smooth = FALSE, color = list(color = "#A5ABB6", highlight = "#FF4F09", arrows = "to")) %>%
  visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE) %>%
  visLayout(randomSeed = 123) %>%
  visInteraction(navigationButtons = TRUE) %>%
  visLegend()


# Customize the background color with custom CSS
network <- network %>% 
  htmlwidgets::onRender("
    function(el, x) {
      el.style.backgroundColor = 'black';
    }
  ")


network

Spinglass Exposing the calibrating parameters for this community detection algo.

spins: The number of spin states. This is a crucial parameter that can indirectly determine the number of communities detected. A larger number of spins allows the algorithm to explore a more significant number of potential community structures.

gamma: The resolution parameter influences the size of the communities detected. Higher values of gamma can lead to finer, smaller communities, while lower values can result in larger, more aggregated communities.

Code

# Choosing the weights
#E(bmi2021_graph)$Value
#E(bmi2021_graph)$betweenness 
#E(bmi2021_graph)$eigenvector 
#E(bmi2021_graph)$closenesss
#E(bmi2021_graph)$degree 

# Assuming bmi2021_graph is already created and is an igraph object
# Calibrating Parameters for Spinglass
spins = 9 # (2 to 30)
gamma = 1 # (1 to 15)
weights = E(bmi2021_graph)$betweenness


spinglass_communities <- cluster_spinglass(bmi2021_graph, 
                                           spins = spins, 
                                           gamma = gamma,
                                           weights = weights)

# Assign community membership to nodes for visualization
V(bmi2021_graph)$community <- as.factor(membership(spinglass_communities))

# Visualize the graph
set.seed(123) # For reproducible layout
g <- ggraph(bmi2021_graph, layout = 'kk') + 
  geom_edge_link(aes(width = Value), alpha = 0.5) +
  scale_edge_width(range = c(0.5, 8)) + 
  geom_node_point(aes(fill = bmi_usd_price, size = bmi_usd_price, color = community), show.legend = FALSE) +
  geom_node_text(aes(label = country), repel = TRUE, size = 2, hjust = -0.2, vjust = 0.2, color = "black") +
  theme_graph() +
  labs(fill = "BMI USD Price", color = "Community", size = "Betweenness")

print(g)

This code performs a sweep over a range of spins and gamma values for the spinglasscommunity detection algorithm, calculates modularity scores for each, and visualizes these scores to identify the optimal spins | gamma value for community detection in your graph. The ggplot2 visualization will help you easily spot the steps value that maximizes the modularity score, indicating an optimal community structure.

Code

library(igraph)

gammas <- seq(1, 30, by = 5)
spins <- seq(2, 30, by = 5)

# Initialize a matrix to store modularity scores
modularity_matrix <- matrix(nrow = length(spins), ncol = length(gammas), dimnames = list(spins, gammas))

# Grid search over gamma and spins
for (g in gammas) {
  for (s in spins) {
    # Spinglass community detection
    spinglass_communities <- cluster_spinglass(bmi2021_graph, spins = s, gamma = g, weights = E(bmi2021_graph)$eigenvector)
    # Calculate modularity
    mod_score <- modularity(spinglass_communities, weights = E(bmi2021_graph)$eigenvector)
    # Store the modularity score
    modularity_matrix[which(spins == s), which(gammas == g)] <- mod_score
  }
}

# Visualization
library(ggplot2)
library(reshape2)

Attaching package: 'reshape2'

The following object is masked from 'package:tidyr':

    smiths

Code

# Melt the matrix for ggplot
modularity_df <- melt(modularity_matrix, varnames = c("Spins", "Gamma"), value.name = "Modularity")

# Heatmap of modularity scores
ggplot(modularity_df, aes(x = Gamma, y = Spins, fill = Modularity)) +
  geom_tile() +
  scale_fill_gradient2(low = "red", high = "blue", mid = "white", midpoint = median(modularity_df$Modularity), space = "Lab", name="Modularity") +
  theme_minimal() +
  labs(x = "Gamma", y = "Spins", title = "Modularity Scores for Different Gamma and Spins in Spinglass Community Detection") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Calling out the communities - These communities will be passed on for regression.

Code

# Assuming you've already created 'node_stats' as previously discussed
community_membership <- membership(spinglass_communities)
node_names <- V(bmi2021_graph)$country  # Ensure this matches your node naming

if(length(community_membership) == length(node_names)) {
    node_stats <- data.frame(
        node = node_names,
        community = as.factor(community_membership)
    )
  
    # Group by community and concatenate country names
    countries_in_community <- aggregate(node ~ community, data = node_stats, FUN = function(x) paste(x, collapse = ", "))
  
    # Print the list of countries by community
    print(countries_in_community)
} else {
    print("Mismatch in lengths of community memberships and node names.")
}

   community                                          node
1          1                                 Brazil, Egypt
2          2                          South Korea, Hungary
3          3                              Spain, Australia
4          4                              Portugal, Turkey
5          5                              Uruguay, Croatia
6          6                                Viet Nam, Oman
7          7                                      Thailand
8          8                        Czech Republic, Jordan
9          9 Canada, Hong Kong, India, Malaysia, Indonesia
10        10                                       Belgium
11        11                   Switzerland, Sweden, Israel
12        12                         United States, Kuwait
13        13          United Kingdom, Qatar, Poland, China
14        14                                   Netherlands
15        15                            Chile, Philippines
16        16                               Singapore, Peru
17        17                                    Costa Rica
18        18                  Denmark, South Africa, Ghana
19        19                               Bahrain, Mexico
20        20                     New Zealand, Saudi Arabia
21        21                 Lebanon, United Arab Emirates
22        22                                        France
23        23                                Japan, Ukraine
24        24                                        Greece
25        25                                       Germany
26        26                             Argentina, Russia
27        27                                       Austria

Code

# Perform Spinglass community detection
communities_spinglass <- cluster_spinglass(bmi2021_graph, weights = E(bmi2021_graph)$Value)

# Add community info to the nodes dataframe
bmi_node_2021$community <- communities_spinglass$membership

# Ensure 'id' column for joining and 'label' for visNetwork node labels
bmi_node_2021$label <- bmi_node_2021$country

bmi_node_2021$community <- as.factor(bmi_node_2021$community)

num_communities <- length(unique(bmi_node_2021$community))

# Generate a color for each community
# Use colorRampPalette to create more colors if needed
color_palette <- colorRampPalette(RColorBrewer::brewer.pal(min(12, num_communities), name = "Set1"))
community_colors <- color_palette(num_communities)

# Map community IDs to colors
color_mapping <- setNames(community_colors, levels(bmi_node_2021$community))
bmi_node_2021$color <- color_mapping[bmi_node_2021$community]

bmi_edges_2021_aggregated <- bmi_edges_2021 %>%
  left_join(bmi_node_2021, by = c("Origin" = "country")) %>%
  rename(from = id) %>%
  left_join(bmi_node_2021, by = c("Partners" = "country")) %>%
  rename(to = id) %>%
   group_by(from, to) %>%
    summarise(Value = n()) %>%
 filter(from!=to) %>%
  #filter(Value > 1) %>%
  ungroup()
bmi_node_2021$label <- bmi_node_2021$country

# Create the network visualization
network <- visNetwork(nodes = bmi_node_2021, edges = bmi_edges_2021_aggregated) %>%
  visNodes(aes(fill = "community", color = "community"), shape = "dot", scaling = list(label = list(enabled = TRUE)),
           font = list(size = 20, face = "arial", color = "white", bold = "true")) %>%
  visEdges(smooth = TRUE, color = list(color = "#A5ABB6", highlight = "#FF4F09", arrows = "to")) %>%
  visOptions(highlightNearest = TRUE, nodesIdSelection = TRUE) %>%
  visLayout(randomSeed = 123) %>%
  visInteraction(navigationButtons = TRUE) %>%
  visLegend()

  # Add custom CSS to change the background color
network <- network %>% 
  htmlwidgets::onRender("
    function(el, x) {
      el.style.backgroundColor = 'black';
    }
  ")

network

1. Walktrap

Insight: Walktrap finds communities by simulating random walks on the network, with the idea that walks are “trapped” within densely connected parts of the graph. Implication for Big Mac Index: The communities detected could reflect clusters of countries where trade flow is more frequent or robust, possibly due to geographic proximity or trade agreements, which in turn may affect the local pricing strategies included in the Big Mac Index.

2. Spinglass

Insight: This method uses a spin model from statistical mechanics and is particularly good at detecting community structures where communities may be hierarchically nested. Implication for Big Mac Index: Spinglass communities may reflect complex and layered trade relationships, possibly uncovering hierarchies in economic interactions that could influence the Big Mac Index.

3. Comparison Summary:

The Walktrap method may reveal more about regional trading blocks or groups of countries with stronger internal trade links that could have similar Big Mac pricing due to shared market conditions. The Spinglass method can uncover more nuanced and multi-level communities, possibly highlighting intricate trade relationships and their cascading effects on the Big Mac Index.

3.4 Model Evaluation

Code

library(igraph)
#E(bmi2021_graph)$Value
#E(bmi2021_graph)$betweenness 
#E(bmi2021_graph)$eigenvector 
#E(bmi2021_graph)$closenesss
#E(bmi2021_graph)$degree 
# Walktrap


weights = bmi_node_2021$Value
steps = 30 # (100 to 500)

walktrap_communities <- cluster_walktrap(bmi2021_graph, 
                                         steps = steps,
                                         weights = weights)

wt_modularity <- modularity(walktrap_communities)
print(paste("Walktrap Modularity:", wt_modularity))

[1] "Walktrap Modularity: 0.606179045040309"

Code

# Spinglass

spins = 9 # (2 to 30)
gamma = 1 # (0 to 10)
weights = E(bmi2021_graph)$betweenness 


spinglass_communities <- cluster_spinglass(bmi2021_graph, 
                                           spins = spins, 
                                           gamma = gamma,
                                           weights = weights)

sg_modularity <- modularity(spinglass_communities)
print(paste("Spinglass Modularity:", sg_modularity))

[1] "Spinglass Modularity: 0.104565961782889"

Code

# Decision based on higher modularity
if (wt_modularity > sg_modularity) {
  print("Walktrap is more optimized based on modularity.")
} else {
  print("Spinglass is more optimized based on modularity.")
}

[1] "Walktrap is more optimized based on modularity."